function [feaScore, feaIdx] = LLEscore(X, K, gamma)
% z: the weight of feature, the samller, the better
% LLE score: a new filter-based unsupervised feature selection method based on nonlinear manifold embedding and its application to image recognition
% Chao Yao, Ya-Feng Liu, Member, IEEE, Bo Jiang, Jungong Han, and Junwei Han, Senior Member, IEEE
% 
% implemented by Liang Du (csliangdu@gmail.com)
%

if ~exist('K', 'var')
    K = 5;
end

if ~exist('gamma', 'var')
    gamma = 1e-5;
end

[nSmp, nFea] = size(X);

M = lle_quad(X, K, gamma);
M = full(M);
z = zeros(nFea,1);
for iDim = 1:nFea
    Mi = lle_quad(X(:, iDim), K, gamma);
    z(iDim) = sum(sum( (M - Mi).^2));
end
feaScore = z;
[~, feaIdx] = sort(feaScore, 'ascend');
end


function M = lle_quad(X, K, gamma)
if ~exist('K', 'var')
    K = 5;
end

if ~exist('gamma', 'var')
    gamma = 1e-5;
end
[nSmp, nDim] = size(X);

X2 = sum(X.^2, 2);
distance = repmat(X2, 1, nSmp) + repmat(X2', nSmp, 1) - 2 * (X * X');
[sorted, index] = sort(distance, 2, 'ascend');
neighborhood = index(:, 2:(1+K));

val = zeros(nSmp, K);
opts = [];
opts.Display = 'iter';
for i1 = 1:nSmp
    xi = X(i1, :);
	Xk = X(neighborhood(i1,:)',:);
    x = quadprog(Xk * Xk' + gamma * eye(K), - Xk * Xi',[],[],ones(1,K),1,[],[],[],opts);
    val(i1, :) = x';
end
W = zeros(nSmp);
W(sub2ind([nSmp, nSmp], repmat([1:nSmp]', K, 1), neighborhood(:))) = val(:);
M = eye(nSmp) - W - W' - (W' * W);
M = (M + M')/2;
end


